
folders;
addpath('../../Toolbox/export_fig/');

%% Set preferences for country grid 
% which cities to keep?
popthreshold = 0.1;       % only load places with population above threshold
pick_largest = 0.6;       % =1 to use most populated place in a cell to locate the node in the cell, otherwise choose centroid
max_N_cities = inf;       % keep the max_N_cities largest cities
N_cat = 10;                % for plotting only. number of population quantiles is N_pop_cat+1.
figures = 0;

% stock all code control variables
code_control.popthreshold = popthreshold;       
code_control.pick_largest = pick_largest;        
code_control.max_N_cities = max_N_cities;    
code_control.N_cat = N_cat;   

% which group to look at?
group = 1; % 1-Mercosur, 2-ANDES

%{ 
Argentina, Brazil, Paraguay and Uruguay (Mercosur)
Bolivia, Peru and Colombia, Ecuador (ANDES)
%}

% Choose 
if group == 1
    country_list = {'AR', 'BR',  'PY', 'UY'}; %, 
    filepaths = {'argentina_adminstrative_province/argentina_adminstrative_province.shp', ...
        'bra_adm_ibge_2020_shp/bra_admbnda_adm1_ibge_2020.shp',  ...
        'pry_adm_dgeec_2020_shp/pry_admbnda_adm1_DGEEC_2020_2.shp', ...
        'ury_adm_2020_shp/ury_admbnda_adm1_2020_2.shp'};

elseif group == 2
    
    country_list = {'PE', 'CO', 'EC', 'BO'}; %, , 'BO'
    filepaths = {'per_adm_ign_20200714_shp/per_admbnda_adm1_ign_20200714.shp', ...
                 'col-administrative-divisions-shapefiles/col_admbnda_adm1_mgn_20200416.shp', ...
                 'ecu_adm_inec_20190724_shp/ecu_admbnda_adm1_inec_20190724.shp', ...
                 'bol_adm_gb2014_shp/bol_admbnda_adm1_gov_2020514.shp'}; %
   
    %{
    country_list = {'CO', 'EC'}; %, , 'BO'
    filepaths = {'col-administrative-divisions-shapefiles/col_admbnda_adm1_mgn_20200416.shp', ...
                 'ecu_adm_inec_20190724_shp/ecu_admbnda_adm1_inec_20190724.shp'}; %
                  %}
end
% determine the connected country folder based on the group of countries
if group == 2
    datafolder_connected_countries =  './connected_countries_data2/'  ;
end
         
Ncountries = length(country_list);

% First use administrative borders to construct grid cells ,
% https://data.humdata.org/dataset/paraguay-administrative-level-0-boundaries?
% Initialize empty structure for polygons
states = struct('X', [], 'Y', [], 'BoundingBox', [], 'Country', [], 'Geometry', 'Polygon');

k = 1;
for country_n=1:length(country_list)
        
    country_icc = char( country_list(country_n) );  % country ICC code
    country = icc2name( country_icc );           % country name

    datafolder_LAC_country =  [ datafolder_LAC,country_icc,'/'];
                
    shp =  shaperead([datafolder_LAC_country,  cell2mat(filepaths(country_n))]);
    
    if strcmp(country_icc, 'PE') % take out Iquitos
        % First find where they do not intersect
        % Then dissolve that area back with the original area 
         [lat,lon] = polysplit( shp(21).Y',shp(21).X' ); % split the country into its polygons

         if length(lat)>1
             polyareas = [];
             for j=1:length(lat)
                 polyareas(j) = polyarea( lon{j},lat{j} );
             end

            % May need to adjust this for Latin American countries with
            % islands?
             [~,keep_poly] = maxk( polyareas,1 );
             [ shp(21).Y ,shp(21).X ] = polyjoin( lat(keep_poly),lon(keep_poly) );
             shp(21).Y = shp(21).Y';
             shp(21).X = shp(21).X';            
         end

        % add new shape back to original shapefile 
    end
    
    if strcmp(country_icc, 'EC') % take out Iquitos
        % First find where they do not intersect
        % Then dissolve that area back with the original area 
        for i=1:length(shp)
         [lat,lon] = polysplit( shp(i).Y',shp(i).X' ); % split the country into its polygons

             if length(lat)>1
                 polyareas = [];
                 for j=1:length(lat)
                     polyareas(j) = polyarea( lon{j},lat{j} );
                 end

                % May need to adjust this for Latin American countries with
                % islands?
                 [~,keep_poly] = maxk( polyareas,1 );
                 [ shp(i).Y,shp(i).X] = polyjoin( lat(keep_poly),lon(keep_poly) );
                 shp(i).Y = shp(i).Y';
                 shp(i).X = shp(i).X';   
             end
        end
        % add new shape back to original shapefile 
    end
    
    if strcmp(country_icc, 'BO') % take out Iquitos
        % First find where they do not intersect
        % Then dissolve that area back with the original area 
         [lat,lon] = polysplit( shp(4).Y',shp(4).X' ); % split the country into its polygons

         if length(lat)>1
             polyareas = [];
             for j=1:length(lat)
                 polyareas(j) = polyarea( lon{j},lat{j} );
             end

            % May need to adjust this for Latin American countries with
            % islands?
             [~,keep_poly] = maxk( polyareas,1 );
             [ shp(4).Y,shp(4).X] = polyjoin( lat(keep_poly),lon(keep_poly) );
             shp(4).Y = shp(4).Y';
             shp(4).X = shp(4).X';   
         end

        % add new shape back to original shapefile 
    end
        
    if strcmp(country_icc, 'AR')
        shp(1, :) = [];
        shp(length(shp), :) = [];        
    end
        
    if strcmp(country_icc, 'EC') % take out the Galapagos Islands 
        shp(9, :) = [];
        shp(7, :) = [];        
    end    
    
    if strcmp(country_icc, 'PE') % take out Iquitos
       shp(16, :) = [];     
    end  
    
    % for Peru and Bolivia, we have to deal with the lake problem
    % using the country borders and the shapefile governing the region in
    % which the lake exists, the gridmap cell will be the union of these
    % areas.    

    if strcmp(country_icc, 'BO') % take out Iquitos
        bo =  load([ datafolder_LAC, 'BO/','country_bounds.mat']);
    end
    
    for row=1:length(shp)
        [ row length(shp) ]

        Y = shp(row).Y;
        X = shp(row).X; 

     %  if ~strcmp(country_icc, 'EC') % take out the Galapagos Islands 
            [Y1, X1, cerr, tol] = reducem(Y', X', 0.01);
            shp(row).Y = Y1';
            shp(row).X = X1'; 
      % end

        states(k).X = shp(row).X;
        states(k).Y = shp(row).Y;
        states(k).BoundingBox = shp(row).BoundingBox;
        states(k).Country = country_icc;
        states(k).Geometry = 'Polygon';       

        k=k+1;
    end    
end

mapshow(states)

% Load in the resulting files 
% load Population Data
file  = [datafolder_connected_countries,'popplaces_sedac.mat'];   % load SEDAC
temp = load(file);   
places_sedac = temp.places_sedac;
clear temp

% total population
totpop = sum( cell2mat( {places_sedac.population} ) );

%%  %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% bring income data %%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%

file  = [ datafolder_connected_countries,'places_gecon.mat'];
load( file,'places_gecon' );

%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% bring altitude data %%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

file  = [datafolder_connected_countries,'relief_etopo.mat'];
load( file,'relief_etopo' );     

%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% discretize the map %%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nquantile=10;
[ places_grid,gridmap,edges,places_gecon ] = connected_country_grid(states ,places_sedac, relief_etopo,places_gecon, Nquantile, group );
         
% From the list of neighbors, generate a list of edges with coordinates
% create map structure
Xtemp = cell2mat({places_grid.X});
Ytemp = cell2mat({places_grid.Y});
unique_edges = edges.unique_edges;
for j=1:length( unique_edges )
    
    % coordinates of each link
   discretized_roads(j).Geometry = 'Line';
   discretized_roads(j).X = [ Xtemp( unique_edges(j,1) ) Xtemp( unique_edges(j,2) ) NaN ];
   discretized_roads(j).Y = [ Ytemp( unique_edges(j,1) ) Ytemp( unique_edges(j,2) ) NaN ];
   
end

% export unique edges to OSM 
coordinates = zeros( length(unique_edges), 4 );  % [ X1 Y1 X2 Y2 ]
for j=1:length( unique_edges )

    % coordinates of each link
    coordinates( j, 1:4 ) = [  Xtemp( unique_edges(j,1) )  Ytemp( unique_edges(j,1) ) ...
                            Xtemp( unique_edges(j,2) )  Ytemp( unique_edges(j,2) ) ];
end

% write coordinates to spreadsheet (one per country)
writetable(array2table(coordinates), [path_raw_mapdata, '/OSMR/connected_coordinates', num2str(group),'.csv']);

%%%% NOTE: TO PROCEED, YOU MUST RUN osm_routing_connected.R 
% TO GENERATE TRAVEL TIMES BETWEEN LOCATIONS for each group 
% executed = system('R CMD BATCH osm_routing_connected.R');

% run file osm_routing.R file to get travel speeds 
speeds = readtable([path_raw_mapdata, '/OSMR/connected_coordinates', num2str(group),'_wtimes.csv']);

for j=1:length( unique_edges )
    
        % coordinates of each link
       if ~isnan(cell2mat({speeds.speed(j)}))
           discretized_roads(j).totdist = cell2mat({speeds.travel_distance(j)});
           discretized_roads(j).avI = cell2mat({speeds.speed(j)});
           discretized_roads(j).keep = 1;

       else
           discretized_roads(j).totdist = deg2km(distance(coordinates(j, 1), coordinates(j, 2), ...
               coordinates(j, 3), coordinates(j, 3)));
           discretized_roads(j).avI = 1e-4;
           discretized_roads(j).keep = 0;  

       end
   
       discretized_avI(j) = discretized_roads(j).avI;
       discretized_dist(j) = discretized_roads(j).totdist;
   
end

%% Given discretized roads, export graph for solver
places2 = places_grid;
graph_export.J = length(places2);

x = zeros(graph_export.J,1);
y = zeros(graph_export.J,1);
for j=1:length(places2)
    
    graph_export.nodes{j}.neighbors = places2(j).neighbors; 
    graph_export.nodes{j}.x = places2(j).X;
    graph_export.nodes{j}.y = places2(j).Y;
    
    x(j) = places2(j).X;
    y(j) = places2(j).Y;
    
end

% vectors with coordinates of each location
graph_export.x=x;
graph_export.y=y;

% graph with discretized network
EdgeTable = table( unique_edges,discretized_avI',discretized_dist', ...
                'VariableNames',{'EndNodes','Weight','Distance'} );   % 'Weight' is the average investment         

grid_graph = graph( EdgeTable );

% adjacency matrix - defined based on avI (variable 'Weight' defined above)
graph_export.adjacency = full( adjacency( grid_graph ) );    

% export some matrixes
nn = numnodes( grid_graph );
[s,t]=findedge( grid_graph );

% matrix of avI
avI_mat = full( sparse( s,t,grid_graph.Edges.Weight,nn,nn ) );  % this gives an upper triangular matrix
avI_mat = avI_mat+avI_mat.' - diag( diag( avI_mat ) );          % this filles the lower-triangular part
graph_export.avI = avI_mat;

% matrix of distances
distance_mat = full( sparse( s,t,grid_graph.Edges.Distance,nn,nn ) );       % this gives an upper triangular matrix
distance_mat = distance_mat+distance_mat.' - diag( diag( distance_mat ) );  % this filles the lower-triangular part
graph_export.distance = distance_mat;

%% add matrix with bilateral distances    

J = graph_export.J;
all_distances = zeros( J,J );
Xtemp = cell2mat({places2.X});
Ytemp = cell2mat({places2.Y});

for i=1:J
        all_distances( i,: ) = deg2km( distance( ones( 1,J )*Ytemp(i),ones( 1,J )*Xtemp(i),...
                                                 Ytemp,Xtemp ) );
end

% distances
graph_export.all_distances = all_distances;
    
% population 
graph_export.L = cell2mat({places2.population})';

% income
for i = 1:length(places2)
    places2(i).income = cast(places2(i).income,'double');
end
    
graph_export.Y = cell2mat({places2.income})';

% geography
graph_export.av_alt = cell2mat({places2.av_alt})';
graph_export.sd_alt = cell2mat({places2.sd_alt})';
graph_export.rugged = cell2mat({places2.rugged})';

% Export data to graph for calibration
country_graph.places_grid=places_grid;
country_graph.gridmap=gridmap;
country_graph.places2=places_grid;
country_graph.discretized_roads = discretized_roads;
country_graph.graph_export = graph_export;
country_graph.edges = edges;  
country_graph.unique_edges = unique_edges;      
country_graph.code_control = code_control;

% Save the grid 
if group ==1
    country='MERCOSUR';
elseif group==2
    country='ANDES';
end

save( [ path_save_grids,country, '_grid.mat' ],'country_graph' );
    